Discontinuities of the exchange-correlation kernel and charge-transfer excitations in 

time-dependent density functional theory 



Maria Hellgren and E. K. U. Gross 
Max-Planck-InsUtut fur Mikrostrukturphysik, Weinberg 2, D-06120 Halle, Germany 

(Dated: February 16, 2012) 

We identify the key property that the exchange-correlation (XC) kernel of time-dependent density 
functional theory must have in order to describe long-range charge-transfer excitations. We show 
that the discontinuity of the XC potential as a function of particle number induces a space -and 
frequency-dependent discontinuity of the XC kernel which diverges as r — ^ oo. In a combined donor- 
acceptor system, the same discontinuity compensates for the vanishing overlap between the acceptor 
and donor orbitals, thereby yielding a finite correction to the Kohn-Sham eigenvalue differences. This 
mechanism is illustrated to first order in the Coulomb interaction. 



I. INTRODUCTION 

The theoretical prediction of the excitation spectra 
of interacting electronic systems is a major challenge 
in quantum chemistry and condensed matter physics. 
A method that has been gaining popularity in the 
past years is time-dependent density functional the- 
ory (TDDFT) offering a rigorous and computationally 
efficient approach for treating excited states of large 
molecules and nanoscale systems. In TDDFT the in- 
teracting electronic density is calculated from a system 
of non-interacting electrons moving in an effective lo- 
cal Kohn-Sham (KS) potential [1]. The KS potential is 
the sum of the external, the Hartree and the so-called 
exchange-correlation (XC) potential, Wxc(rt), virhich, due 
to the Runge-Gross theorem [2], is a functional of the 
density. In the linear response regime, the excitation en- 
ergies can be extracted from the poles of the linear den- 
sity response function. As a consequence, given the vari- 
ational derivative /xc(i"t, r't') = Svy^d'^t) / 6n{Y't')^ also 
known as the XC kernel [3], it is possible to formulate 
an RPA-like equation for the exact excitation spectrum 
[4, 5]. In practical calculations, this equation is solved 
using some approximate XC potential and kernel where 
the most popular ones are based on the adiabatic local 
density approximation (ALDA), leading to kernels local 
in both space and time. Despite this simple structure op- 
tical excitations of small molecules are successfully pre- 
dicted. However, several shortcomings have also been re- 
ported: cxcitons in solids arc not captured, [6, 7] double 
excitations are missing [8] and charge-transfer (CT) exci- 
tations are qualitatively incorrect [9-12]. In this work we 
will be concerned with the last problem and to see why 
ALDA fails in this case, we consider a charge transfer be- 
tween two neutral Coulombic fragments. The asymptotic 
limit of the excitation energy is then given by 

UCT =Id- Aa - 1/R, (1) 

where Id is the ionization energy of the donor, Aa is the 
affinity of the acceptor and R is their separation. In 
TDDFT the starting point is the KS system which yields 
the exact Id but only an approximate Aa- Thus, the XC 
kernel must both account for the 1/R correction and shift 



the KS affinity. The linear response equations involve, 
however, only matrix elements of /xc between so-called 
excitation functions $iQ(r) = (pi{r)ipa{r), i.e., products 
of occupied and unoccupied KS orbitals. As the distance 
between the fragments increases these products vanish 
exponentially and thus there is no correction to the KS 
eigenvalue differences unless the kernel diverges [13, 14]. 
Kernels from the ALDA, or adiabatic GGA's for that 
matter, do not contain such divergency and it is as yet 
not understood how this extreme behavior should be in- 
corporated in approximate functionals. 

Whenever two subsystems are spatially well-separated 
it is possible to treat one of the subsystems in terms of an 
ensemble containing states with different number of par- 
ticles. DFT has been generalized to non-integer charges 
and as an important consequence it was found that the 
XC potential jumps discontinuously by a constant at in- 
teger particle numbers in order to align the highest oc- 
cupied KS eigenvalue with the chemical potential [15]. 
In this way, the true affinity. A, is equal to the sum of 
the KS affinity. As, and the discontinuity. Not surpris- 
ingly, it has therefore been argued that the discontinuity 
must play an important role in describing CT excitations 
within TDDFT [16]. 

So far, only the XC potential has been the target for in- 
vestigating discontinuities in DFT and TDDFT [17-19]. 
In this work we instead examine possible discontinuities 
of the XC kernel. We demonstrate the existence of a 
discontinuity and we study its properties. Furthermore, 
we give an explicit example through a numerical study 
of the EXX functional. Finally, as a first application, 
we demonstrate the crucial role of the discontinuity for 
capturing CT excitations in linear response TDDFT. 



II. DERIVATIVE DISCONTINUITY IN DFT 

We start by considering a static system of electrons 
described by a statistical operator p = J2k'^k\'^k){'^k\, 
where |^'fc) denotes the ground state of k particles corre- 
sponding to the Hamiltonian H = T + V + J dr w{r)n(r), 
in which T is the kinetic energy, V the inter-particle inter- 
action and w is the external potential. The ground-state 
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energy Eq of the system with average number of parti- 
cles N is obtained by minimizing the functional [n] = 
F[n]+ J drw{r)n{r), where F[n] = minp^„ Tr[p{f + V)], 
under the constraint that N = Jdrn(r). At the mini- 
mum n = n[w,N] coincides with the ground-state den- 
sity. The XC energy is defined as i?xcM = — ?s[7i] — 
U[n\, where Ts[n] = minp_s.„ Tr[/5T] is the non-interacting 
kinetic energy and U[n\ is the Hartree energy. Assum- 
ing that the density can be reproduced by an ensemble of 
non-interacting electrons the XC part of the KS potential 
is given by Wxc(r) ~ (5i?xc/<5n(r). In general Eq and, in 
particular, i?xc['^[w, A^]] has derivative discontinuities at 
integer particle numbers A^o [20, 21]. The partial deriva- 
tive of Exc with respect to A, 



dE^ 
dN 



dr iJxc(r) 



dn{r) 
dN 



(2) 



has two sources of discontinuous behavior: (i) the quan- 
tity 



_ dn{r) 

^ dN ' 



(3) 



known as the Fukui function, may have different right 
and left limits, /+ and (the superscript ± refers to the 
value of the quantity at A^ = Aq + O^), (ii) the XC poten- 
tial may be discontinuous, v^^{r) = v^^{r) + Axe, where 
Axe is a constant. Below we will show that also the sec- 
ond variation of E^c with respect to the density, the static 
XC kernel /xc(r, r') = 6vxc{'r)/Sn{r'), has discontinuities 
which are related to derivative discontinuities in the den- 
sity itself. As pointed out in previous work [5, 22, 23], 
the particle number conserving density response is unaf- 
fected by adding to /xc a function depending only on one 
of the coordinates. In line with the results of Ref. 24 we 
therefore argue that the discontinuities of /xc must be of 
the form 



/+ (r, r') = /-,(r, r') + .9xc(r) + gxc(r'). 



(4) 



In the following we will show a simple procedure for de- 
termining Axe and (7xc(r) which is useful whenever i?xc 
is an implicit functional of the density via, e. g., the KS 
Green function. 



A. XC potential 

For N > No, we write fxc(i") = ^xd''^) + ^xc(r) and 
cast Eq. (2) into 



/.rAx.(r)/(r)^^-/ 



- drt--,(r)/(r). 



(5) 



In the limit A^ Nq , Axc(r) — > Axe and we find a 
formal expression for the discontinuity of v^c 



A. 



dEx 



dN 



drvxcir)f^ir). 



(6) 



This expression can be used as the starting point for de- 
riving the well-known MBPT-formula for the correction 
to the gap [25], as we will now demonstrate. From the 
Klein functional within MBPT [26] it is possible to con- 
struct an XC energy functional in terms of the KS Green 
function G's(r,r',a;) [27]. In this case Sxe = ^Exc/SGs, 
where Exe is the self-energy evaluated at Gs ■ The deriva- 
tive of i?xc with respect to A^ is then given by 



dExc 
dN 



= 1^1 rfrdr'Sxe(r,r',c.) 



dGs{r,r',ij) 
dN 



(7) 



In order to evaluate the derivative of Gs with respect 
to A^ we consider an ensemble described by a spin- 
compensated mixture of states with electron number 



A^o and Aq - 
Ae [AcA^o 



1. The KS ensemble Green function for 
- 1] is given by 



No 



G.(r,r',..) = J2 



V3fe(r)(y9fe(r') 



k=l 

P<^L(r)<^L(r') 



k=No+2 



2 w 



1 



£k + ir] 
p\ (^L(r)<^L(r') 



(8) 



where p ^ N — Nq and the subscript L signifies the lowest 
unoccupied molecular orbital (LUMO) of the KS system, 
which is considered partially occupied and partially un- 
occupied. Notice that the KS orbitals and eigenvalues 
Efc also depend on A^ via the KS potential Vs. The deriva- 
tive of Gs with respect to A^ is now easily carried out 



dGs{r,r',uj) ^ 1 ^i^{r)^i^{r') 1 ^i^{r)^i^{r') 
2 u! — Eh — iv 2 oj — Eh + if] 
6Gs{r,r',io) dVsir,) 



dN 



dri- 



SVsivi) 

From Eq. (9) and Eq. (7) wc find 
dE, 



dN 



(9) 



= J drdrVJ(r)S+(r,r',eL)¥'J(r') 
Jdrdr'dr,^Ur,r',u.^ ^G.(r,r',.) ay.(rO 



(5T4(ri) dN 



(10) 



The second term on the right hand side of Eq. (6) can 
be written as 

J drv-,{r)f+{r) = J dr v-,{r)\ip+{r)\' 

Sn{r') dVsivi) 



dr'dri Wxe(i"') 



SVs (ri) dN 



(11) 



The discontinuity Axe is now easily determined. The 
second terms on the right hand side of Eqs. (10-11) will 
cancel by virtue of the linearized Sham-Schliiter equation 
[28] and we find 



Axe = J rfrdrVL(r)S+(r,r',eL)<y9L(r') 



^^I■Wxe(r)l'PL(r)|^ 



(12) 
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where we have omitted the superscript on the orbitals where we have used the chain rule 
since they are continuous with respect to N. Eq. (12) 
agrees with the one in Rcf. 25. 



B. XC kernel 



Svy,c{r) _ f ^^, Sv^c{r) 5n{r') ^^^^ 



Sw{ri) J Sn{r') 6w{ri) 



Next, we turn to the XC kernel for which we exhibit 
the discontinuities by taking the functional derivative of 
dE^c/dN (Eq. (2)) with respect to w. This yields 



5w{yi) dN 



drdr'x(ri,r')/xc(r',r)/(r) 
6fir) 



5w{ri) 



(13) 



and identified the linear density response function 
x(ri,r') — Sn{r')/Sw{ri). Then, we write /xc(r,r') = 
/^(r, r') + 5xc(i", I"')' insert in Eq. (13), and take the 
limit N Nq . According to the discussion above, in 
this limit (7xc(r, r') — > 5xc(r) + 5xc(r')- Using further- 
more that / c?r x(r, r') = and / drf{r) = 1 we find the 
following equation for the discontinuity gxc of /xc: 



J drx(ri,r)gxc(r) 



S dE^ 



5w{vi) dN 



J drdr'x(ri,r)/x-(r,r')/+(r') - J drv+{r) 



(15) 



We notice that Eq. (15) only determines ^xc up to a constant. This constant can, however, easily be fixed by 
considering d^E^c/dN'^ in the limit N —J- Nq 



(16) 



2 J rfr/+(r)gxc(r) = ^ ^- J dr' v+{r') ^ ^- J drdr' f+{r)f-{ry)f+{r') 



This equation does not allow an arbitrary constant in gxc- 
Consequently, Eq. (15) and Eq. (16) together uniquely 
determines the discontinuity of /xc. 

To gain insight on the r-dependence of 5xc(i") we 
employ a common energy denominator approximation 
(CEDA) [29] to Eq. (15). To do so, we first notice that 
the derivative with respect to w can be replaced by the 
derivative with respect to Vs using the chain-rule since 
the density is a functional of w only via Vs . The CEDA 
then allows us to partially invert the KS response func- 
tion Xs analytically. We will focus on the left hand side 
of Eq. (15) and on the last term on the right hand side. 
These terms are less sensitive to the approximation used 
for E'xc and should therefore give rise to a general behav- 
ior. If all energy denominators are set to the constant Ae 
we find on the left hand side of Eq. (15) 



J rfi"Xs(ri,i-).gxc(r) ~ -^n(ri),9xc(ri) 

+ J^y"*7(ri,r)5xc(r)7(r,ri) (17) 

where 7 is the KS density matrix. If we focus on the last 
term in Eq. (15) and use /"''(r) ~ |93L(r)P we find in the 



CEDA 



Sf{r) 



-^kL(ri)|2wxc(i"i) 



6Vsiri 



drv+{r) 



+ ^</'L(ri) J dr7(ri,r)(pL(r)wxc(r) (18) 

From Eq. (17) and Eq. (18) we can extract an approxi- 
mate asymptotic behavior 



= (r) 



\Mr)\ 
n(r) 



= (r) 



(19) 



where / is the ionization energy [30] and Ag the KS affin- 
ity. Thus, we can conclude that if / > g^c contains a 
term which diverges exponentially as r — 00. That the 
first term of Eq. (15) would exactly cancel the term of 
Eq. (19) is highly unlikely and we will explicitly see that 
within an approximation that accounts for the derivative 
discontinuity as the EXX approximation this is not the 
case. To obtain Eq. (19) we have used the CEDA but we 
will show below that the discontinuity obtained from the 
full solution of Eq. (15) exhibits the same behavior. In 
addition, we will demonstrate that this feature is respon- 
sible for sharp peak structures as well as divergences in 
the kernel of a combined donor-acceptor system. 
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III. DISCONTINUITY OF THE DYNAMICAL 
XC KERNEL 

So far the analysis has been Hmited to the static case. 
To investigate if the kernel has discontinuities at finite 
frequency the discussion above must be generalized to 
an ensemble which allows the number of particles to 
change in time. To this end, we consider the follow- 
ing statistical operator p{t) = '^k'^k{t)\'^ k[t)){'^ k{t)\^ 
where ak{t) are given time-dependent coefficients whose 
sum is equal to 1 and Y^kit)) is the many-body state of 
k particles at time t. It is possible to prove a Runge- 
Gross-like theorem for this ensemble which allows us 
to define the XC potential as a functional of the en- 
semble density. In this way, the functional derivative 
5vy^c{'^t) / 5n{Y't') contains no arbitrariness, but leaves the 
possibility for a discontinuity of the form /^(r, r'; t—t') = 
/- (r, r'; t - t') + g^v; t-t')+ g^v'; t-t'). If we vary 
'yxc(rt) with respect to the time-dependent number of 
particles N{t') and evaluate the derivative at the ground 
state with N = Nq an expression for the frequency- 
dependent discontinuity (7xc(r, i^) can be derived 



(5wxc(ri) 



5N{t') 



= J dr'/xc(r,r',t-t')/+(r') 



+ j dr'gxc(r', t - t')J+{v') + g^v, t - t'), (20) 

where we again have used the chain rule. Below we will 
see with a numerical example how the frequency depen- 
dence modifies the discontinuity making the divergency 
stronger than in the static case and allows for a cor- 
rect description of CT excitations in a combined donor- 
acceptor system. This will be illustrated in the time- 
dependent (TD) EXX approximation. For a derivation 
and analysis of the EXX kernel we refer the reader to Ref. 
[31, 32]. From now on the subscript xc will be replaced 
by X to denote quantities in the TDEXX approximation. 



IV. RESULTS 



In the single-pole (SP) approximation [4] of TDDFT 
the XC correction to the KS excitation energy u!q 
is given by twice the matrix element {q\fxci'-o)\q) — 
J (ir(ir'$g(r)/xc(r, r', i^)$<j(r') at w = uig, where the in- 
dex q = ia corresponds to an arbitrary excitation. In Ref. 
[33] it has been shown that in TDEXX 2{q\f^{ijjq)\q) = 
(a[I]x — Wx|a} — (*|Sx — Wxl*} — {aa\v\ii), where v is the 
Coulomb interaction. Considering a CT excitation be- 
tween HOMO (i = H) and LUMO (a = L) we see 
that the last term goes as 1/R. Setting, as usual [34], 
(H|I]x — Ux|H) = and using the previously mentioned 
result (L|Sx — WxlL) = Ax, we can deduce that [35] 

wcT = WHL +2(HL|w + /x(a;HL)|HL) 

-> tJHL + Ax - l/R. (21) 




5 10 15 20 25 15 20 25 
X [a.u.] X [a.u.] 

FIG. 1: Left: EXX potential and kernel for a heteronuclear 
system of four electrons. Right: the discontinuity g,c{x, 0) cal- 
culated from Eq. (15) for the isolated two-electron subsystem 
as well as G{x) for A'^ = 2.0001. Note that the potentials have 
been rescaled and shifted for better visibility. 



The kernel /x thus produces a finite correction if evalu- 
ated at WHL and yields exactly the results corresponding 
to first order Gorling-Levy perturbation theory [33, 36- 
38]. In the following we will see that it is the discon- 
tinuity of the kernel that yields this correct result. We 
have deliberately used the SPA and not the full solution 
of Casida equations in conjunction with the EXX kernel, 
a procedure which would imply the inclusion of higher 
orders in the explicit dependence on the Coulomb inter- 
action. Our motives are to study an exact property of 
the kernel well captured in the SPA of TDEXX but may 
be subject to errors inherent to the approximation when 
including higher orders. 

We model [39] a stretched diatomic molecule in terms 
of ID atoms described by Q/y^ {x — xq)^ + 1, where xq 
is the location of the atom and Q is the nuclear charge, 
and replace everywhere the Coulomb interaction v with 

+ T. We study 



r'\2 



a soft-Coulomb interaction l/-\/(a 
two different systems, one ionic and one neutral system. 
In the ionic system the discontinuity is important already 
at the level of the XC potential whereas in the neutral 
system the discontinuity appears only in the XC kernel. 



A. Ionic system 

In the first example we study an ionic system and set 
Q = 2 on the left atom (donor) and that of the right 
atom (acceptor) to Q = 4 and solve the ground-state KS 
problem with 4 electrons. In the ground state at inter- 
nuclear separation i? = 10 a.u. we find two electrons 
on each atom and in Fig. 1 (left panel) the EXX po- 
tential is displayed (fade line) in arbitrary units. Two 
steps are clearly visible, one between the atoms and an- 
other one on the right side of the acceptor. As a con- 
sequence, v-x is shifted upwards in the acceptor region 
placing the KS LUMO of the isolated acceptor above 
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the HOMO of the isolated donor. This imphes that the 
KS afhnity of the acceptor becomes closer to the true 
affinity. In the same figure and panel we also display 
the quantity, Fiii^{x,uj) = / (ia;'/x(a;, a;', w)3>HL(a;'), at 
w = 0. The function i^jjL is seen to have peaks in corre- 
spondence with the steps of and is shifted downwards 
over the acceptor with respect to the donor. Despite the 
fact that $HL tends to zero as R increases the peaks of 
-Fhl become sharper and higher, and the shift increases 
in size. The right panel in the same figure shows the 
function G{x) = / dx' f^{x, x' ,0)f{x'), accessible from 
Eq. (15), for the isolated acceptor when TV = 2.0001 
(full line), as well as the discontinuity g^{r,0) in the 
limit = 2+ (dashed line). The potentials are also 
shown (fade lines) calculated from the same ensembles, 
i.e., N = 2.0001 and 2. The function G{x) has 

peaks whose positions follow the steps of Ux and whose 
hight increases as N approaches 2+. In the same limit, 
also the difference between the G(oo) value and the value 
of G in the central plateau-like region increases consis- 
tently with the fact that / dx G{x)f{x) has to remain 
finite (see Eq. (16)). Eventually G{x) turns into the 
discontinuity g-xi^, 0) which diverges as x — > oo, in agree- 
ment with our previous analysis. Notice that CEDA has 
not been made here. If we compare Fhl and G from 
the different panels we see a very similar structure. We 
therefore conclude that the peaked structure of the kernel 
in the donor-acceptor system is just the discontinuity of 
the ensemble EXX kernel. As discussed above, most part 
of the CT excitation energy is contained already in the 
KS eigenvalue differences due to the step in iix- A more 
critical example is therefore the neutral system, studied 
below. 



B. Neutral system 



T ' I ' I ' ' I ' r 




To ' 20 ' 30 ' 10 20 30 

X [a.u.] X [a.u.] 



FIG. 2: Same system as in Fig. 1 but with six electrons. Left: 
AEXX kernel. Right: TDEXX kernel at cjhl. Note that the 
potentials have been rescaled for better visibility. 
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FIG. 3: CT excitation energies as a function of separation R 
in different approximations. 



An example where the kernel needs to account for Ax 
is the same system but with 6 electrons, i.e., a neutral 
system. The ground-state has 2 electrons on the left atom 
(acceptor) and 4 electrons on the right atom (donor). In 
Fig. 2 we plot the EXX potential (fade line) for i? = 12 
a.u. and a step-like structure between the atoms can be 
observed. However, as R is increased the step reduces 
in size and eventually goes to zero. In the left panel we 
plot Fhl at w = 0, i.e., in the adiabatic (AEXX) approx- 
imation, and for different separations R. Again, we find 
a peak structure in the kernel between the donor and 
the acceptor as well as a shift that increases exponen- 
tially with R. In this case, compared to the previous, the 
peaks are less pronounced but the step due to the plateau 
is much larger. Evaluating the kernel at the first KS CT 
excitation frequency increases the exponential growth of 
the step by around a factor of two (right panel). Thus, 
whereas the overall shape remains unaltered the magni- 
tude of the shift is strongly influenced by the frequency 
dependence. This fact plays a crucial role in the descrip- 
tion of the CT excitation for this system. We notice here 



that even if the step in Wxc disappears in the dissociation 
limit the step in /xc remains. This is not a contradiction 
as the discontinuity might show up only in the second 
derivative. Fig. 3 illustrates the behavior of the SP CT 
excitation energy as a function of R for the system of Fig. 
2 in four different approximations. In TDEXX with the 
correction (HL|/x(whl)|HL) we find that the divergency 
of the kernel over the acceptor exactly compensates for 
the decreasing overlap $hl, thus yielding a finite value as 
i? — > oo as well as the right 1 /R asymptotic behavior as 
it should according to Eq. (21). In the adiabatic case we 
find instead that (HL|/x(0)|HL) tends to as i? — s- cxd, al- 
though reproducing the fully frequency dependent result 
up to around R = 8. We notice that even if the kernel 
is very large over the acceptor it will not affect the exci- 
tations which are localized there since any constant will 
vanish by the fact that integrates to zero. Thus only 
excitations involving a transfer of charge from one atom 
to the other will be influenced by the discontinuity. 
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V. CONCLUSIONS 

In conclusion we have analyzed the discontinuity of the 
XC kernel of an ensemble with time-dependent particle 
numbers. In a combined system of two atoms we have 
seen that the divergency of the discontinuity as r cxo 



can generate a kernel which diverges in the dissociation 
limit, and thus compensate for the vanishing overlap of 
acceptor and donor orbitals. This feature is crucial for 
the description of CT excitations but may also be im- 
portant whenever there are excitations for which the KS 
orbital overlap is too small to give a correction. 
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